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We have computed the one-body density matrix pi in solid ^He at T = K using the Shadow 
Wave Function (SWF) variational technique. The accuracy of the SWF has been tested with an 
exact projector method. We find that off-diagonal long range order is present in pi for a perfect hep 
and bcc solid *He for a range of densities above the melting one, at least up to 54 bars. This is the 
first microscopic indication that Bose Einstein Condensation (BEC) is present in perfect solid *He. 
At melting the condensate fraction in the hep solid is 5 x 10~^ and it decreases by increasing the 
density. The key process giving rise to BEC is the formation of vacancy-interstitial pairs. We also 
present values for Leggett's upper bound on the superfiuid fraction deduced from the exact local 
density. 

PACS numbers: 67.80.-s 



The supersolid state, a solid with superfiuid proper- 
ties, is moving out from theoretical speculations as a re- 
sult of the observation of non classical rotational inertia 
(NCRI) in solid ''He in Vycori and very recently in the 
bulk-. The initial theoretical suggestion^ii of the super- 
solid state was based on the possible presence of vacancies 
in the ground state of a Bose quantum solid. In addi- 
tion these vacancies have to be mobile in order to give 
rise to Bose Einstein Condensation (BEC) of ''He atoms. 
Experimentally no evidence has been found for the pres- 
ence of vacancies at very low temperature and this is in 
agreement with the results of microscopic theory^'^ which 
gives an energetic cost of about 15 K for the formation of 
a vacancy in bulk solid *He. However it was almost im- 
mediately recognized^ that the presence of ground state 
vacancies is only one possible mechanism for NCRI and 
what is really needed is that atoms are not localized at 
the lattice sites but are delocalized via exchange or other 
processes. This gives the possibility of having in the wave 
function (wf) a phase which governs collectively the mo- 
tion of the atoms. The existence of a supersolid state in 
^He is therefore strictly related to the question of local- 
ization or delocalization of particles and, of course, this is 
a topic of general interest. This is the case, for instance, 
of cold alkali atoms in a periodic potential* . Experiments 
with ''He give access to the superfiuid fraction Ps/ p which 
turns out to be^ at most of order 2%. Theoretically only 
upper bounds on p^ have been obtained up to nowSiifliii 
and no microscopic theory has given evidence for a su- 
persolid phase. The commonly accepted view is that one 
can have a finite ps in a three dimensional system if there 
is BEC in the system so a central quantity to compute is 
the off diagonal one-body density matrix pi{r, r '), whose 
Fourier transform represents the momentum distribution. 

In this article we address the computation of pi for 
sohd ^Hc at T = K based on a variational wf, a shadow 
wf (SWF). In a previous computationiS we have found 
that the presence of vacancies in the solid induces a BEC 
which is proportional to the concentration of vacancies. 
On the other hand the large energy of formation of a va- 
cancy makes the probability of having such defects at low 



temperature vanishing small. Here we study in the per- 
fect solid the large distance behavior of pi(r, r'), specif- 
ically if pi has a non zero limit at large distance (off 
diagonal long range order, ODLRO) which implies BEC. 
With perfect solid we mean that the number of maxima 
in the local density p{r) is equal to the number of ^He 
atoms. The use of SWF is specially useful in the present 
context because with such wf the crystalline order is an 
effect of the spontaneously broken symmetry so that lo- 
cal disorder processes like exchange of two or more parti- 
cles, creation of vacancy-interstitial pairs (VIP) or more 
complex processes are in principle allowed. The major 
finding of our computation is the presence of a small but 
finite condensate for a range of densities above melting. 
The variational theory is very useful to describe strongly 
interacting systems like liquid or solid ^He but it is al- 
ways open to debate how much the results depend on the 
ansatz on the wf, specially for quantities other than the 
energy. In order to give indication on the reliability of 
our SWF we present some results on quantities like the 
degree of local order, of localization and of the local den- 
sity obtained also from an exact computation basedi^ on 
a projection algorithm (SPIGS), a Path Integral Ground 
State'^'' method which uses a SWF as the starting wf. 

In a SWF the correlations between atoms are intro- 
duced both explicitly by a Jastrow factor and also in 
an implicit way by coupling with a set of subsidiary vari- 
ables, which are called "shadow" variables'^5 (one shadow 
for each quantum particle), which are integrated over. 
All expectation values are computed by a Monte Carlo 
(MC) method and the statistical sampling of maps 
the quantum system of N particles in a system of N spe- 
cial interacting triatomic "molecules'— which consist of 
a ^He atom and two shadows. The accuracy of the SWF 
technique is well documented and it has been possible to 
treat also disorder phenomena in a quantum solid, like a 
vacancj^j^ or even the interfacial region between a solid 
and a liquid at coexistencei^. As functional form for the 
correlating factors contained in the SWF we have taken 
the ones used in ref^; as interatomic interaction we have 
used a standard Aziz potential. 
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TABLE L Upper bounds for the superfluid fraction in 
hep bulk soUd ''He eomputed at different densities with the 
SWF technique, the SPIGS method and the GM for p{r). P 
is the pressure from the SPIGS equation of state, a is the 
standard deviation of the Gaussians used in the GM. 
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A SWF can be interpreted as a first projection step in 
imaginary time of a Jastrow wf via a variationally op- 
timized imaginary time propagator^^. With SPIGS one 
goes beyond the variational theory by adding successive 
projection steps in the imaginary time propagation with 
the full Hamiltonian and in this way we are able to com- 
pute exact expectation values on the ground state with- 
out extrapolations. With these two Quantum MC meth- 
ods no a priori equilibrium positions for the solid phase 
are required, the Bose symmetry is manifestly maintained 
and atoms can be delocalized. 

The equation of state given by SWF is in good agree- 
ment with the results of SPIGS, for instance at melting 
(p = 0.029 A-3) the energy per particle is -5.12 K, 0.76 
K above the SPIGS resullii^. Also for a vacancy there is 
an excellent agreement, the formation energy at melting 
is 15.7 K with SPIGS and 15.6 K with SWF^ and the 
jumping rate is similar in the two computations. In ad- 
dition to the energy one would like to know the accuracy 
of SWF in describing the microscopic local processes and 
to this end we have computed the local density p{r) and 
the static structure factor S{k). We find again agreement 
between SWF and SPIGS with SWF giving a slightly 
more ordered state. For instance at melting the main 
Bragg peak of the hep solid is about 17% higher than 
the SPIGS result. Detailed comparison will be presented 
elsewhere and here we focus only on the results for the up- 
per bound for the superfluid fraction Ps/ p obtained by 
Leggetlji. This bound depends on the averaged density, 
p[z) = J d^p{r), where z is a longitudinal coordinate and 
^ is a suitable set of transversal coordinates^. We have 
chosen as z axis the one which gives the lowest value of 
f^; in the hep crystal this is the FA direction which is 
perpendicular to the basal plane. It should be noticed 
that usually in Quantum MC computations the center 
of mass of the system is not fixed and this would alter 
the local density specially around the minima. Therefore 
when we compute p{'r) we have modified the sampling al- 
gorithm to keep the center of mass of the system fixed. In 
Tab^we show the upper bounds, f^, for the superfluid 
fraction obtained for hep bulk solid ^He with the SPIGS 
and SWF methods. There is substantial agreement, the 



variational being always lower as a consequence of the 
larger degree of local order. A popular representation of 
p(r) is the one in terms of sum of Gaussians centered on 
the lattice sites. We have fitted our p{f) with this Gaus- 
sian Model (GM) by using the standard deviation a as 
fitting parameter. We find that the GM gives an excel- 
lent representation of the integrated density along planes, 
what is needed in the Leggett's inequality, the deviation 
being below 4%. The resulting bounds given by GM with 
a fitted on the SPIGS p{r) are also shown in TabU The 
bound given by SWF is similar to the value computed 
previousljiS with the GM fitted on a different variational 
theory. In RefiiSiii a lower /+ has been obtained by us- 
ing a better variational ansatz: the phase of the wf is a 
function of r and not only of the longitudinal coordinate 
z as in Ref.^. This bound computed with the GM for a 
given crystal lattice is simply a function of the "localiza- 
tion parameter" Saslow's computation is for an fee 
crystal but if we neglect the difference between the hep 
and the fee lattice, using the a in TabU and the results 
in Ref^, we can estimate an which goes from about 
0.2 at p = 0.029A-3 to about 0.005 at p = 0.044A-3 . 
A SPIGS computation for fee at p = 0.029A~'^ gives a a 
in the GM which is only 2% lower than in hep crystal. 
These values of /^^ are compatible with the experiments 
but is about one order of magnitude larger than the 
experimental value of Ps/p and, in any case, it is only 
an upper bound so it is not very conclusive. A word of 
caution on the GM is in order. If this model gives an ex- 
cellent representation for the integrated density p{z), the 
accuracy is lost when we consider the local density p(r): 
in the region of the minima of p(r) deviations greater 
than 100% are found. 

The one-body density matrix pi(r, r') is given by 
the overlap between the normalized many-body ground 
state wf ^{R) and \E'(i?') where configuration R' = 
{r',r2, ■■t'^n} differs from R = {r,r2, ■■j^'n} only by the 
position of one of the atoms in the system; if ^'(-R) is 
translationally invariant as in our case when the center 
of mass is not fixed, pi only depends on the difference 
r — r': 

Pi[r-f') = Njdr2-dfN'i*[R)'^{R') . (1) 

It is possible to interpret the integrand in Eq.(^| as a 
probability density^^; then pi can be computed by sam- 
pling the integrand in Eq.(^ and by histogramming the 
occurrence of the distance d — f — r' . In the following 
we will call "half" particles the particles with coordi- 
nates r and r ' because they have just ^ the correlation 
strength with the other iV — 1 particles (with coordinates 
{r2, rw}) and no direct correlation between them. The 
method of computation has been described in Ref.^S,. Ab- 
sence of ODLRO corresponds to the two "half" particles 
forming a "molecule" whereas presence of ODLRO corre- 
sponds to a finite probability of dissociation up to infinite 
distances. 

We have computed pi{r — r') along the nearest neigh- 
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FIG. 1: pi{f — f') at difTerent densities for hep and bcc 
perfect ''He crystals and for the same crystals with a finite 
concentration Xv of vacancies. 



bors (nn) direction which is FK in a hep and [111] in a 
bcc crystal. In Fig^we report pi for a perfect hep crys- 
tal at different densities at and above melting and for a 
bcc crystal at p = 0.02898A~'^. It is clear that at melt- 
ing and at p = O.OSlA^'^ pi reaches a plateau at large 
distance whereas at p = O.OSSA"'^ pi steadily decreases 
with increasing distance. By averaging the tail in pi for 
distances greater than 14 A, we find at the melting den- 
sity a condensate fraction ric = (5.0 ± 1.7) x 10^^ in a 
perfect hep and (7.6 ± 1.7) x 10^^ in a perfect bcc crys- 
tal. This is the first microscopic indication that BEC is 
present in perfect solid "^Hei^. At p = O.OSlA^'^ we find 
Uc = (2.0 ±0.4) X 10"'^ and at p = O.OBBA'^ the tail is so 
much depressed that the size of the simulation box is too 
small to conclude if ODLRO is present; in this case we 
can only say that the condensate fraction, if any, is lower 
than 10~®. The simulation box which has been used to 
compute pi is cubic and contains N = 432 '^He atoms for 
the bcc crystal, it is elonged in the FK direction and con- 
tains TV = 360 ''He atoms for the hep crystal. By chang- 
ing the size of the boxes at fixed density we have checked 
that our results for pi have no finte size effect within the 
statistical errors of our computations. This is shown in 
the inset of Fig^for the hep crystal at p = 0.029A~'^; 
by averaging the tail in pi for distances greater than 14 
A, we find Uc = (3.9 ± 1.7) x 10"^ for N = 288 and 
Uc = (5.7 ± 2.0) X 10"^ for N ^ 432, to be compared 
with the value given above: Uc = (5.0 ±1.7) x 10~^ for 
N — 360. In refii^ we found that a finite concentration 
of vacancies induces a condensate fraction which de- 
pends linearly on Xv In Fig^we show also pi when a 
vacancy is present both in hcpi^ and bcc. Taking into 
account the value of x^ of the computation we estimate 
that at p = 0.029A~'^ the condensate fraction due to a 
finite concentration of vacancies is equal to the one in a 
perfect hep crystal when x^ — 1.5 x 10"^. 




FIG. 2: Projection of 100 successive configurations of the 
particles (white circles) and the two "half" particles (black 
circles) in a basal plane of an hep crystal at p = 0.03lA~"^. 
pi is the probability distribution of the two "half" particles. 



All the computed pi show oscillations which reflect the 
crystalline order in the system. However these oscilla- 
tions are not the same in the perfect and in the defected 
solid. We find that when a vacancy is present the maxima 
of the oscillations in pi correspond to multiples of the nn 
distance c?nn- This is an indication that in presence of 
a vacancy the main mechanism which contributes to the 
separation of the two "half" particles is that one of them 
moves through the crystal following the vacancy which 
is very mobilc^*^. The different positions of the maxima 
in pi for the perfect crystal suggest a different micro- 
scopic process for the ODLRO in this case. By analysing 
the particle configurations sampled in our runs we find 
that the secondary peak of pi, located at about 5 A, 
always corresponds to a configuration in which the two 
"half" particles occupy two nn lattice positions slightly 
distorted by the presence of one interstitial ''He atom be- 
tween them. In Fig|21we show 100 successive configura- 
tions of the particles which corresponds to this event. In 
the formalism of second quantization pi (r — r ') is equal to 
the expectation value of the composite event where one 
^He atom is destroyed at f ' and one is created at r; then 
it is possible to interpret the event in Fig|21 as the cre- 
ation of a VIP. The same process is found in bcc crystal. 
After this first step the two "half" particles have a finite 
probability of moving away one from another by exchange 
processes with the other atoms and this gives rise to the 
other maxima of pi at larger distance. By analysing the 
particle configurations corresponding to these other max- 
ima we find that a VIP is present in all the configurations. 
Similar processes were considered in Ref.^^ as a necessary 
condition for the supersolid phase, but there it is argued 
that VIP cannot be present. Our results disagree with 
this hypothesis. In order to characterize the anisotropy 
of pi as function of d we have computed pi when the two 
"half" particles are no more constrained to lie on the nn 
direction but can freely move in a plane. In Fig|31one can 
see that pi in a perfect bcc crystal is strongly anisotropic 
for distances up to about 6 A, and the maxima of pi are 
in the direction of nn. However at greater distances pi 
becomes nearly isotropic and we conclude that our esti- 
mation of the BEC fraction is not affected by the pre- 
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FIG. 3: pi(f-f ') at p = 0.02898A~^ for a perfect bcc crystal 
with f — r' lying in the plane [-101]. 

vious restriction on d. Similar results are obtained for 
hep. It is interesting to notice that when a vacancy is 
present the anisotropy of pi persists up to greater dis- 
tances (data not shown). Also in FigQlone can see that 
the oscillations of pi are more persistent with increasing 
distance in the crystal with a vacancy; this is another 
indication that different microscopic processes are at the 
origin of the ODLRO in the perfect and in the defected 
solid ^He. The exchange of atoms and VIPs are present 
not only in pi but also in j^"!^. At melting about ev- 
ery 2x10'^ MC steps an atom has a displacement larger 
than dnn and in many cases this is associated with the 
presence of an interstitial. In principle one can devise an 



algorithm based on SPIGS to compute pi exactly. How- 
ever at present this appears to be a major computational 
problem. In any case we have given solid evidence that 
SWF overestimates the degree of local order so that we 
should expect that the SWF results for the BEG fraction 
are an underestimation of the exact values. 

In conclusion we have shown that solid *He at T = 
K has BEG at melting density and above at least up to 
54 bars whereas we find a vanishing BEG at 90 bars. 
Thus BEG should be at the basis of the NGRI observed 
experimentally^. Our result has been obtained from an 
advanced variational theory the accuracy of which has 
been tested with a projector method on the exact ground 
state. The key process giving rise to ODLRO is the for- 
mation of a VIP. Such defects have a finite probability 
to be present in the ground state of the system; they are 
not permanent excitations but simply rare fluctuations of 
the perfect crystal induced by the large zero-point mo- 
tion. In other words the number of atoms is equal to the 
number of lattice sites and, at the same time, atoms are 
delocalized. Since the ground state is the vacuum of the 
elementary excitations of the system we conjecture that 
a branch of low energy excitations different from phonons 
should be present in solid ''He. Such excitations should 
have an important role in determining the critical tem- 
perature. It is a possibility that this branch is related to 
some experimental results which have been interpreted 
in term of an excitation with energy of about 2 KSSi. 
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